# Librarieslibrary(tidyverse)library(gmRi)library(rnaturalearth)library(scales)library(sf)library(gt)library(patchwork)library(lmerTest)library(emmeans)library(merTools)library(tidyquant)library(ggeffects)library(performance)library(gtsummary)library(gt)library(sizeSpectra)library(ggdist)library(pander)conflicted::conflict_prefer("select", "dplyr")conflicted::conflict_prefer("filter", "dplyr")conflicted::conflicts_prefer(lmerTest::lmer)# Processing functionssource(here::here("Code/R/Processing_Functions.R"))# ggplot themetheme_set(theme_gmri(rect =element_rect(fill ="white", color =NA), strip.text.y =element_text(angle =0),axis.text.x =element_text(size =8),legend.position ="bottom"))# vectors for factor levelsarea_levels <-c("GoM", "GB", "SNE", "MAB")area_levels_long <-c("Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight")area_levels_all <-c("Northeast Shelf", area_levels)area_levels_long_all <-c("Northeast Shelf", area_levels_long)# table to join for swapping shorthand for long-hand namesarea_df <-data.frame(area =c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "All"),survey_area =c("SS", "GoM", "GB", "SNE", "MAB", "Northeast Shelf"),area_titles =c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "Northeast Shelf"))# Degree symboldeg_c <-"\u00b0C"# Function to get the Predictions# Drop effect fits that are non-significant ###get_preds_and_trendsignif <-function(mod_x){ modx_preds <-as.data.frame(# Model predictionsggpredict( mod_x, terms =~ est_year * survey_area * season) ) %>%mutate(survey_area =factor(group, levels = area_levels),season =factor(facet, levels =c("Spring", "Fall")))# Just survey area and year modx_emtrend <-emtrends(object = mod_x,specs =~ survey_area*season,var ="est_year") %>%as_tibble() %>%mutate(zero =0,non_zero =if_else(between(zero, lower.CL, upper.CL), F, T))# Preds with signif modx_preds %>%left_join(select(modx_emtrend, survey_area, season, non_zero))}
Storycrafting - Northeast Shelf Spectra
I think we need to step back a bit and find the main story again–which is about how size structure of the NES fish community has changed over time, whether that varies by region, and what covariates are associated with those changes.
The Aug draft you developed really focused on comparing regression models to consider drivers based on the fish community size spectrum as viewed from the perspective of the Wigley species set vs. the full species set.
However, I didn’t see figures of those spectra in that draft (maybe there are somewhere else and I’ve lost track). The recent slide deck really focuses on seasonal differences in the size spectra and influence of temperature. - KMills
Starting from Square Uno
I’d like for us to start by rebuilding the story of how the size structure has changed over time and regional differences.
As such, I would like for us to look at plots of the size spectrum slope for (1) the all-species length spectrum and (2) the Wigley species biomass spectrum from annual, seasonal and regional perspectives.
I am thinking this would be two panel plots (all-spp length, Wigley-spp biomass) with 5 subpanels each (NES, GOM, GB, SNE, MAB) that are set up to show lines by seasonal and annual time steps, much like the plot on the right of slide 7 in your most recent slide deck…except with NES added and an annual line added. Keeping the significant trend lines in is helpful. This is the first step, so when you have that ready, please share it with the group. - KMills
Code
# Data used for Wigley estimateswigley_in <-read_csv(here::here("Data/processed/wigley_species_trawl_data.csv"))finfish_in <-read_csv(here::here("Data/processed/finfish_trawl_data.csv"))
Aside: Biomass Fraction in Each Community
If we need to make a decision about which community we want to use, here is the different biomass totals for each:
Code
# glimpse(wigley_in)# Load raw# General tidying only, removal of strata outside our study area, Spring and Fall only# (removes inshore and rarely/inconsistently sampled strata etc.)# shrimps removedtrawl_basic <-read_csv(here::here("Data/processed/trawl_clean_all_data.csv"))# Biomass totalsraw_totals <- trawl_basic %>%distinct(est_towdate, cruise6, station, stratum, tow, comname, biomass_kg) %>%pull(biomass_kg) %>%sum()wig_biomass_total <- wigley_in %>%distinct(est_towdate, cruise6, station, stratum, tow, comname, biomass_kg) %>%pull(biomass_kg) %>%sum()ffish_biomass_total <- finfish_in %>%distinct(est_towdate, cruise6, station, stratum, tow, comname, biomass_kg) %>%pull(biomass_kg) %>%sum()# Make a table for displaytibble("Community"=c("All (including crustaceans)", "All finfish", "Wigley Species Subset"),"Biomass Total"=c(raw_totals, ffish_biomass_total, wig_biomass_total)) %>%mutate(`% of Total`= (`Biomass Total`/ raw_totals)*100) %>%gt()
# Left side:# All species species length spectrum # color = season vs. annual# facet_wrap(~survey_area)# Right side:# Wigley species biomass spectrum# color = season vs. annual# facet_wrap(~survey_area)# Plot over observed data# Contrast seasonal differences(lenspectra_trend_p | bmspectra_trend_p) +plot_layout(guides ="collect")
Median Length & Weight Trends
I also think it would be valuable to build out the story of size change further, so am thinking that plots like those of average (or median) length (all species) and weight (wigley species) for NES and the sub-regions, like those on p. 12 of the attached file above (the very first draft from late 2022) would be useful. - KMills
Code
# Load the median weight/length datawigley_size_df <-read_csv(here::here("Data/processed/wigley_species_size_summary.csv"))finfish_size_df <-read_csv(here::here("Data/processed/finfish_species_length_summary.csv"))# shelf-wide summaryfinfish_lengths_shelf <-read_csv(here::here("Data/processed/shelfwide_finfish_species_length_summary.csv"))wigley_sizes_shelf <-read_csv(here::here("Data/processed/shelfwide_wigley_species_size_summary.csv"))# Join themsize_results <-bind_rows(list("Finfish Community"=bind_rows(finfish_size_df, finfish_lengths_shelf),"Wigley Subset"=bind_rows(wigley_size_df, wigley_sizes_shelf)), .id ="community") %>%mutate(survey_area =factor(survey_area, levels = area_levels_all) )# Get trends:len_mod_ffish <-lmer(formula = med_len_cm ~ est_year * survey_area * season + (1|est_year),data = finfish_size_df)# len_mod_wig <- lmer(# formula = med_len_cm ~ est_year * survey_area * season + (1|est_year),# data = wigley_size_df)wt_mod_wig <-lmer(formula = med_wt_kg ~ est_year * survey_area * season + (1|est_year),data =drop_na(wigley_size_df, med_wt_kg))# Get the predictions and flag whether they are significantsize_trend_predictions <-bind_rows(list("med_len_cm-Finfish Community"=get_preds_and_trendsignif(len_mod_ffish),#"med_len_cm-Wigley Subset" = get_preds_and_trendsignif(len_mod_wig),"med_wt_kg-Wigley Subset"=get_preds_and_trendsignif(wt_mod_wig) ), .id ="model_id") %>%separate(model_id, into =c("metric", "community"), sep ="-") %>%mutate(metric =fct_rev(metric))# Left side: # Median length - all species# color = season vs. annual# facet_wrap(~survey_area)size_long <- size_results %>%pivot_longer(cols =c(med_len_cm, med_wt_kg), names_to ="metric", values_to ="val") %>%mutate(survey_area =fct_relevel(survey_area, area_levels_all),season =fct_rev(season))# just length for scaleslen_plot <- size_long %>%filter(metric =="med_len_cm", community =="Finfish Community") %>%drop_na(val) %>%ggplot() +geom_ribbon(data =filter(size_trend_predictions, metric =="med_len_cm", non_zero ==TRUE),aes(x, ymin = conf.low, ymax = conf.high, fill = season),linewidth =1, alpha =0.35) +geom_point(#data = filter(size_long, metric == "med_len_cm"),aes(est_year, val, color = season),size =0.8, alpha =0.5) +geom_ma(aes(est_year, val, color = season, linetype ="5-Year Moving Average"), n =5, ma_fun = SMA) +geom_line(data =filter(size_trend_predictions, metric =="med_len_cm", non_zero ==TRUE),aes(x, predicted, color = season, linetype ="Significant Trend"),linewidth =1, alpha =1) +scale_fill_gmri() +scale_color_gmri() +facet_grid(survey_area ~ metric*community,scales ="free") +theme(strip.text.y =element_blank()) +labs(x ="Year",y ="Length (cm)",color ="Season",linetype ="Linetype",fill ="Season")# Right: # Median weight - wigley species# color = season vs. annual# facet_wrap(~survey_area)# weight plotswt_plot <- size_long %>%filter(metric =="med_wt_kg") %>%drop_na(val) %>%ggplot() +geom_ribbon(data =filter(size_trend_predictions, metric =="med_wt_kg", non_zero ==TRUE),aes(x, ymin = conf.low, ymax = conf.high, fill = season),linewidth =1, alpha =0.35) +geom_point(aes(est_year, val, color = season),size =0.8, alpha =0.5) +geom_ma(aes(est_year, val, color = season, linetype ="5-Year Moving Average"), n =5, ma_fun = SMA) +geom_line(data =filter(size_trend_predictions, metric =="med_wt_kg", non_zero ==TRUE),aes(x, predicted, color = season),linewidth =1, alpha =1) +scale_fill_gmri() +scale_color_gmri() +facet_grid(survey_area ~ metric*community,scales ="free") +labs(x ="Year",y ="Weight (kg)",color ="Season",linetype ="Linetype",fill ="Season")(len_plot | wt_plot) +plot_layout(guides ="collect")
There are a number of years when median weight surges in MAB. For these situations I think we likely could trace these down to exceptionally high catches of larger species like elasmobranchs.
We also don’t see a shelf-wide decline in size so I might need to check what Friedland did or prepare for a confrontation.
Small/Large Fish Index
The following figures frame “large fish” as individuals with bodymass greater than the 90th percentile for all years within that region. “Small fish” are considered smaller than the 10th percentile of the mass distribution.
The following table shows what those percentiles are for each region.
NOTE: Percentiles are determined using DescTools::Quantile() which uses numlen_adj to weight everything by abundance.
Code
# # Get 95th percentile as threshold# DescTools::Quantile(# wigley_in$length_cm, # weights = wigley_in$numlen_adj, # probs = c(0.05, 0.1, 0.5, 0.90, .95)) %>% # t() %>% # as_tibble()# Do them separatelyarea_size_quants <-bind_rows(mutate(wigley_in, survey_area ="Northeast Shelf"), wigley_in) %>%split(.$survey_area) %>%map_dfr(function(x,y){ DescTools::Quantile(# x$length_cm, x$ind_weight_kg,weights = x$numlen_adj, probs =c(0.05, 0.1, 0.5, 0.90, .95)) %>%t() %>%as_tibble() }, .id ="survey_area") %>%mutate(survey_area =factor(survey_area, levels = area_levels_all)) %>%arrange(survey_area)# Show the tablearea_size_quants %>%mutate_if(is.numeric, round, 2) %>%gt() %>%tab_header(title ="Wigley Species - Body Mass Quantiles (kg)")
# Plot them as loing-term average# Plot the SFIsfi_p <- lfi %>%ggplot() +geom_point(aes(est_year, SFI, color = season),size =0.8, alpha =0.5) +geom_ma(aes(est_year, SFI, color = season, linetype ="5-Year Moving Average"), n =5, ma_fun = SMA) +scale_fill_gmri() +scale_color_gmri() +scale_y_continuous(labels =label_percent()) +theme(strip.text.y =element_blank()) +facet_grid(survey_area~., scales ="free") +labs(x ="Year",y ="Percent of Total Abundance",title ="Small Fish Index\n(< 10th Percentile Bodymass)",fill ="Season",color ="Season",linetype ="")# Plot the LFIlfi_p <-ggplot(lfi) +geom_point(aes(est_year, LFI, color = season),size =0.8, alpha =0.5) +geom_ma(aes(est_year, LFI, color = season, linetype ="5-Year Moving Average"), n =5, ma_fun = SMA) +scale_fill_gmri() +scale_color_gmri() +scale_y_continuous(labels =label_percent()) +facet_grid(survey_area~., scales ="free") +labs(x ="Year",y ="",title ="Large Fish Index\n(> 90th Percentile Bodymass)",fill ="Season",color ="Season",linetype ="")# Plot together(sfi_p | lfi_p) +plot_layout(guides ="collect")
Code
# Plot them as loing-term average# Plot the SFIsfi_p <- lfi %>%ggplot() +geom_point(aes(est_year, small_abund, color = season),size =0.8, alpha =0.5) +geom_ma(aes(est_year, small_abund, color = season, linetype ="5-Year Moving Average"), n =5, ma_fun = SMA) +scale_fill_gmri() +scale_color_gmri() +scale_y_continuous(labels =label_number()) +theme(strip.text.y =element_blank()) +facet_grid(survey_area~., scales ="free") +labs(x ="Year",y ="Total Abundance",title ="Small Fish Index\n(< 10th Percentile Bodymass)",fill ="Season",color ="Season",linetype ="")# Plot the LFIlfi_p <-ggplot(lfi) +geom_point(aes(est_year, large_abund, color = season),size =0.8, alpha =0.5) +geom_ma(aes(est_year, large_abund, color = season, linetype ="5-Year Moving Average"), n =5, ma_fun = SMA) +scale_fill_gmri() +scale_color_gmri() +scale_y_continuous(labels =label_number()) +facet_grid(survey_area~., scales ="free") +labs(x ="Year",y ="",title ="Large Fish Index\n(> 90th Percentile Bodymass)",fill ="Season",color ="Season",linetype ="")(sfi_p | lfi_p) +plot_layout(guides ="collect")
Summary of Trends
Let’s start with these pieces, and we can then discuss doing something with the extremes of the size distribution–changes in small vs. large individuals, where we pre-define size bins as you’ve done in slides 7 and 8, or perhaps using an approach similar to Lora’s approach of defining the upper/lower percentiles of size. - KMillz
Part 2: Relationships to External Factors
The following plots and summaries will relate the Wigley species spectra to Temperature and Live Landings
Code
# # Model Dataset for Wigley Species Subsetwigley_bmspectra_model_df <-read_csv(here::here("Data/model_ready/wigley_community_bmspectra_mod.csv"))# Use one for plotswigley_bmspectra_model_df <- wigley_bmspectra_model_df %>%mutate(survey_area =case_when( survey_area =="GoM"~"Gulf of Maine", survey_area =="GB"~"Georges Bank", survey_area =="SNE"~"Southern New England", survey_area =="MAB"~"Mid-Atlantic Bight"),survey_area =factor(survey_area, levels = area_levels_long),season =fct_rev(season))
Bottom Temperature Changes
Bottom temperatures have increased in all four regions. The distribution of seasonal bottom temperatures are such that Gulf of Maine experiences a much smaller difference in seasonal temperature swings than the other regions.
Code
temp1 <- wigley_bmspectra_model_df %>%ggplot(aes(est_year, bot_temp, color = season)) +geom_point(size =0.8, alpha =0.5) +geom_ma(size =1, n =5, ma_fun = SMA, linetype =1) +scale_color_gmri() +facet_wrap(~survey_area, ncol =1) +scale_y_continuous(labels =label_number(suffix ="\u00b0C")) +labs(y ="Bottom Temperature", x ="Year", color ="5-Year Smooth")# Plot the distribution as a raincloud plottemp2 <- wigley_bmspectra_model_df %>%ggplot(aes(fct_rev(survey_area), bot_temp, color = season, fill = season)) + ggdist::stat_halfeye(alpha =0.4,adjust = .5, width = .6, .width =0, justification =-.2, point_colour =NA ) +geom_boxplot(fill ="transparent",width = .15, outlier.shape =NA ) +## add dot plots from {ggdist} package ggdist::stat_dots(## orientation to the leftside ="left", size =1,## move geom to the leftjustification =1.12, ## adjust grouping (binning) of observations binwidth = .25) +scale_color_gmri() +scale_fill_gmri() +scale_y_continuous(labels =label_number(suffix ="\u00b0C")) +coord_flip() +theme(legend.position ="none") +labs(y ="Bottom Temperature", x ="", color ="Season")temp1 | temp2
The overlap in bottom temperature distributions is very similar to the overlap in mass spectra exponents.
Code
wigley_b_dist_plot <- wigley_bmspectra_model_df %>%ggplot(aes(fct_rev(survey_area), b, color = season, fill = season)) + ggdist::stat_halfeye(alpha =0.4,adjust = .5, width = .6, .width =0, justification =-.2, point_colour =NA) +geom_boxplot(fill ="transparent",width = .15, outlier.shape =NA) +## add dot plots from {ggdist} package ggdist::stat_dots(## orientation to the leftside ="left", size =1,## move geom to the leftjustification =1.12, ## adjust grouping (binning) of observations binwidth = .005) +scale_color_gmri() +scale_fill_gmri() +coord_flip() +labs(y ="ISD Exponent", x ="",title ="Individual Body Mass Distribution (1g - Max)")# Plot next to bottom temperatureswigley_b_dist_plot | temp2
Which makes me lean towards thinking that the size distribution is responsive to the ambient water temperatures, and perhaps less stationary and responsive to long-term trends in some area.
The following model has been used to determine the relationship between landings and seasonal bottom temperature on b.
lmer(
b ~ survey_area*season*scale(roll_temp) + survey_area * scale(log(total_weight_lb)) + (1|est_year),
data = wtb_model_df)
Code
#### Model 2: Temperature & Fishing ##### Drop NA'swtb_model_df <-drop_na(wigley_bmspectra_model_df, total_weight_lb, bot_temp, b) %>%rename(area = survey_area) %>%left_join(area_df) %>%# Do rolling BT within a season * regiongroup_by(survey_area, season) %>%mutate(roll_temp = zoo::rollapply(bot_temp, 5, mean, na.rm = T, align ="right", fill =NA)) %>%ungroup() %>%mutate(yr_num =as.numeric(est_year),yr_fac =factor(est_year),survey_area =factor(survey_area, levels = area_levels),season =factor(season, levels =c("Spring", "Fall")),landings = total_weight_lb,yr_seas =str_c(season, est_year))# Make the modelmod2 <-lmer( b ~ survey_area*season*scale(roll_temp) + survey_area*scale(log(total_weight_lb)) + (1|est_year), data = wtb_model_df)# summary(mod2)# check_model(mod2)
Drive Model Sumamry
Code
mod2 %>% gtsummary::tbl_regression() %>%add_global_p(keep = T) %>%bold_p(t =0.05) %>%bold_labels() %>%italicize_levels() %>%as_gt() %>%tab_header(title ="Wigley Community - Bodymass Distribution & Drivers")
Wigley Community - Bodymass Distribution & Drivers
Characteristic
Beta
95% CI
1
p-value
survey_area
0.030
GoM
—
—
GB
-0.13
-0.29, 0.02
0.088
SNE
0.07
-0.06, 0.19
0.3
MAB
0.05
-0.05, 0.15
0.3
season
0.2
Spring
—
—
Fall
0.06
-0.04, 0.17
0.2
scale(roll_temp)
-0.01
-0.13, 0.10
0.8
scale(log(total_weight_lb))
0.08
0.00, 0.15
0.048
survey_area * season
<0.001
GB * Fall
0.17
0.00, 0.35
0.051
SNE * Fall
-0.12
-0.29, 0.06
0.2
MAB * Fall
-0.22
-0.39, -0.06
0.009
survey_area * scale(roll_temp)
0.007
GB * scale(roll_temp)
-0.24
-0.40, -0.08
0.004
SNE * scale(roll_temp)
0.02
-0.14, 0.18
0.8
MAB * scale(roll_temp)
-0.10
-0.27, 0.07
0.3
season * scale(roll_temp)
0.7
Fall * scale(roll_temp)
0.03
-0.11, 0.16
0.7
survey_area * scale(log(total_weight_lb))
0.4
GB * scale(log(total_weight_lb))
-0.06
-0.17, 0.06
0.3
SNE * scale(log(total_weight_lb))
-0.05
-0.17, 0.07
0.4
MAB * scale(log(total_weight_lb))
-0.07
-0.14, 0.01
0.094
survey_area * season * scale(roll_temp)
0.031
GB * Fall * scale(roll_temp)
0.16
-0.03, 0.36
0.10
SNE * Fall * scale(roll_temp)
-0.14
-0.33, 0.06
0.2
MAB * Fall * scale(roll_temp)
0.06
-0.13, 0.26
0.5
1
CI = Confidence Interval
Temperature Relationship
Code
# Clean up the plot:# Plot the predictions over datamod2_preds <-as.data.frame(ggpredict( mod2, terms =~ roll_temp*survey_area*season)) %>%mutate(survey_area =factor(group, levels = area_levels),season =factor(facet, levels =c("Spring", "Fall")))#### Trend Posthoc ####trend_df <-emtrends(object = mod2,~survey_area * season,var ="roll_temp")# Drop effect fits that are non-significant ###mod2_emtrend <-emtrends(object = mod2,specs =~ survey_area*season,var ="roll_temp") %>%as_tibble() %>%mutate(zero =0,non_zero =if_else(between(zero, lower.CL, upper.CL), F, T))# Limit temperature plotting rangeactual_range <- wtb_model_df %>%group_by(season, survey_area) %>%summarise(min_temp =min(bot_temp)-2,max_temp =max(bot_temp)+2,.groups ="drop")# Plot over observed datamod2_preds %>%left_join(actual_range) %>%filter((x < min_temp) == F, (x > max_temp) == F) %>%left_join(select(mod2_emtrend, survey_area, season, non_zero)) %>%filter(non_zero) %>%ggplot() +geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = survey_area), alpha =0.1) +geom_line(aes(x, predicted, color = season), linewidth =1) +geom_point(data = wtb_model_df,aes(roll_temp, b, color = season),alpha =0.4,size =1) +facet_grid(survey_area~., scales ="free") +scale_color_gmri() +scale_x_continuous(labels =label_number(suffix = deg_c)) +labs(y ="Exponent of ISD",title ="Wigley Species, Body Mass Spectra",x ="5-Year Rolling Average Bottom Temperature")
Landings Relationship
Code
# Clean up the plot:# Plot the predictions over datamod2_preds <-as.data.frame(ggpredict( mod2, terms =~ total_weight_lb * survey_area * season)) %>%mutate(survey_area =factor(group, levels = area_levels),season =factor(facet, levels =c("Spring", "Fall")))#### Trend Posthoc ####trend_df <-emtrends(object = mod2,~survey_area * season,var ="total_weight_lb")#summary(trend_df, infer= c(T,T))# Drop effect fits that are non-significant ###mod2_emtrend <-emtrends(object = mod2,specs =~ survey_area*season,var ="log(total_weight_lb)") %>%as_tibble() %>%mutate(zero =0,non_zero =if_else(between(zero, lower.CL, upper.CL), F, T),group =str_c(survey_area, "X", season))# Plot over observed datamod2_preds %>%left_join(select(mod2_emtrend, survey_area, season, non_zero)) %>%filter(non_zero) %>%ggplot() +geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha =0.1) +geom_line(aes(x, predicted, color = season), linewidth =1) +geom_point(data = wtb_model_df,aes(total_weight_lb, b, color = season),alpha =0.4,size =1) +facet_grid(survey_area~., scales ="free") +scale_color_gmri() +scale_x_continuous(transform ="log10", labels =label_log(base =10)) +labs(y ="Body Mass Spectra Slope (b)",title ="Wigley Species, Body Mass Spectra",x ="Total Landings (lb.)")
Side Stories
Everything below is digging into different nuances of the approach or the community composition. These things are of secondary importance and will be moved out.
Spiny Dogfish Sensitivity
Spiny dogfish constitute a large majority fraction of biomass in the MAB during the Spring survey. Their numbers have also been increasing over time. This acts to shallow the spectra as they are one of the larger species routinely sampled by the trawl survey.
If we remove them from the estimation and re-run the slopes this is what changes.
How much can one species shift the spectra, and is this a definitive/conclusive proof that its the dogfish shallowing the slope in MAB?
Code
# # Run MAB with and without spiny dogfish# # Run the year, season, area fits for the filtered data# no_dogs_bmspectra <- group_binspecies_spectra(# ss_input = filter(wigley_in, comname != "spiny dogfish"),# grouping_vars = c("est_year", "season", "survey_area"),# abundance_vals = "numlen_adj",# length_vals = "length_cm",# use_weight = TRUE,# isd_xmin = 1,# global_min = TRUE,# isd_xmax = NULL,# global_max = FALSE,# bin_width = 1,# vdiff = 2)# # Save those# write_csv(# no_dogs_bmspectra,# here::here("Data/processed/wigley_species_nodogfish_bmspectra.csv"))
Code
# Read them in, do some plottingno_dogs_bmspectra <-read_csv( here::here("Data/processed/wigley_species_nodogfish_bmspectra.csv")) %>%mutate(survey_area =factor(survey_area, levels = area_levels))# Format a littleno_dogs_bmspectra <- no_dogs_bmspectra %>%mutate(est_year =as.numeric(est_year))# Join them togetherdfish_results <-bind_rows(list("Wigley Full"= wigley_bmspectra,"Wigley w/o Spiny Dogfish"= no_dogs_bmspectra), .id ="community") %>%mutate(metric ="bodymass_spectra") %>%mutate(survey_area =fct_relevel(survey_area, area_levels))# Get trends for no dogfishnodfish_mod <-lmer(formula = b ~ est_year * survey_area * season + (1|est_year),data = no_dogs_bmspectra)bmspectra_mod_wig <-lmer(formula = b ~ est_year * survey_area * season + (1|est_year),data = wigley_bmspectra)# Put predictionsdogfish_sens_predictions <-bind_rows(list("bodymass_spectra-Wigley Full"=get_preds_and_trendsignif(bmspectra_mod_wig),"bodymass_spectra-Wigley w/o Spiny Dogfish"=get_preds_and_trendsignif(nodfish_mod)), .id ="model_id") %>%separate(model_id, into =c("metric", "community"), sep ="-") %>%mutate(metric =fct_rev(metric))# Plot the differencesggplot() +geom_ribbon(data =filter(dogfish_sens_predictions, non_zero ==TRUE),aes(x, ymin = conf.low, ymax = conf.high, fill = season),linewidth =1, alpha =0.35) +geom_point(data = dfish_results,aes(est_year, b, color = season),size =0.8, alpha =0.5) +geom_line(data =filter(dogfish_sens_predictions, non_zero ==TRUE),aes(x, predicted, color = season),linewidth =1, alpha =1) +scale_fill_gmri() +scale_color_gmri() +facet_grid(survey_area ~ metric*community, scales ="free") +labs(x ="Year",y ="Exponent of Size Spectra",title ="Shift in Spectra Trend(s) with Exclusion of Spiny Dogfish")
Story Thoughts
Below is a table of some papers operating in the same area or around the same topic that are particularly relevant:
Code
shelf_papers <-tribble(~"Author", ~"Year", ~"Conjecture","Friedland et al.", 2024, "Body size has declined despite more abundance & biomass. Competition is preventing larger sizes.","Krumsick", 2020, "Empirical slopes from trawl survey in labrador sea were steeper than theoretical from nitrogen stable isotope ratios")shelf_papers %>%gt()
Author
Year
Conjecture
Friedland et al.
2024
Body size has declined despite more abundance & biomass. Competition is preventing larger sizes.
Krumsick
2020
Empirical slopes from trawl survey in labrador sea were steeper than theoretical from nitrogen stable isotope ratios
There are a couple core ideas about what changes have occurred or have been ocurring over the Northeast Shelf with relevance to the community that is sampled by the trawl survey: